Glucose Metabolism and Oscillatory Behavior of Pancreatic Islets 
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A variety of oscillations are observed in pancreatic islets. We establish a model, incorporating two 
oscillatory systems of different time scales: One is the well-known bursting model in pancreatic j3- 
cells and the other is the glucose-insulin feedback model which considers direct and indirect feedback 
of secreted insulin. These two are coupled to interact with each other in the combined model, and 
two basic assumptions are made on the basis of biological observations: The conductance gn(ATP) for 
the ATP-dependent potassium current is a decreasing function of the glucose concentration whereas 
the insulin secretion rate is given by a function of the intracellular calcium concentration. Obtained 
via extensive numerical simulations are complex oscillations including clusters of bursts, slow and 
fast calcium oscillations, and so on. We also consider how the intracellular glucose concentration 
depends upon the extracellular glucose concentration, and examine the inhibitory effects of insulin. 

PACS numbers: 87.19.Nn, 87.19.-j, 05.45.-a, 05.45.Xt 



I. INTRODUCTION 

Insulin secretion from pancreatic (3 cells is critical for 
glucose homeostasis in blood. Perturbations from the 
basal glucose concentration induce various oscillations 
with different periods in pancreatic (3 cells: fast oscil- 
lations with periods of several seconds QJ, slow oscil- 
lations with periods of several minutes |2(, and ultra- 
dian oscillations with periods of a few hours Q. Each 
type of oscillations has been studied separately, within 
the respective model according to each time scale. For 
example, slow oscillations can be caused by a feedback 
mechanism having slow reaction time . Independently 
of this, a microscopic model for the bursting mecha- 
nism of the membrane potential in a /3 cell was devised 
and investigated while the bursting activity was also 
demonstrated in a simple two-dimensional map 0. In 
particular, the microscopic bursting model, which was 
constructed in a similar manner to the Hodgkin-Huxlcy 
model for nerve excitation [f| , was improved by the intro- 
duction of the potassium ion K + conductance depending 
dynamically on the ATP concentration However, the 
correlations between the bursting mechanism and insulin 
secretion have not been addressed in an appropriate way 
and there still lacks a model describing the complex phe- 
nomena such as clusters of bursting action potentials , 
fast and slow Ca 2+ oscillations § , which arise from com- 
bined behaviors of different time scales. The purpose of 
this study is to propose a combined model, explaining 
both fast and slow oscillations of action potentials. 

On the other hand, recent experimental investigations 
in vivo and in vitro about the oscillations of insulin secre- 
tion have revealed a more complex picture on the molec- 
ular basis. Stimulation of insulin secretion by glucose 
involves a rise in the cytoplasmic concentration of cal- 
cium ions (Ca 2+ ) in (3 cells. This rise essentially re- 



sults from the following sequence of events, as sketched 
in Fig. □ Glucose transported by GLUT-1 and GLUT- 
2 transporters in a (3 cell raises the ratio of adenosine 
triphosphate (ATP) to adenosine diphosphate (ADP), 
which promotes closure of ATP-sensitive K + channels. 
This generates membrane depolarization, urging voltage- 
dependent calcium channels to open. The subsequent 
increase in free cytosolic Ca 2+ then stimulates insulin se- 
cretion . Note also that there are several experimental 
evidences for insulin to inhibit bursting of (3 cells [Tol [TT| : 
in particular, there is a report that insulin activates ATP- 
sensitive K + channels in (3 cells 12]. Based on these ex- 
periments, we assume two possible inhibitory pathways 
of insulin: indirect effects on gating of K(ATP) channels 
by inhibiting glucose transport through GLUT-2 trans- 
porters by insulin and direct effects of insulin to activate 
the channel. Both negative feedback processes of insulin 
to activate the channels may result in oscillatory behav- 
ior. Oscillations of the membrane potential |13| drive 
oscillations of Ca 2+ JljK leading in turn to oscillations of 
insulin secretion [l{|[l(j- In fact, the measured portal- 
vein insulin secretion rate shows periodic behavior of in- 
sulin secretion with successive secretion and rest. The 
frequency of such pulsatile insulin secretion in vivo and 
that in vitro in an isolated perfused pancreas have been 
reported to vary from four to fifteen minutes and from 
six to ten minutes, respectively 0, 0] . 

These experimental results make it desirable to estab- 
lish a model that can describe the whole feedback pro- 
cess on the microscopic level, with the aforementioned 
macroscopic features incorporated. As an attempt to- 
ward such a goal, we propose a model, connecting the 
macroscopic description of glucose regulation 2] and the 
microscopic mechanism of bursting behavior of (3 cells Q 
via the prescription based on experimental observations: 
The conductance gx(ATP) for the K(ATP) current, i.e., 
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Direct Pathway 



FIG. 1: Mechanism of insulin secretion. Activation and in- 
hibition of GLUT-1 and GLUT-2 transporters by secreted 
insulin are represented by the solid (+) and dashed ( — ) ar- 
rows. Thick arrows describe physical transport of materials 
(glucose and ions). 



the K + current through the ATP-dependent K + channel, 
is taken to be a decreasing function of the intracellular 
glucose concentration whereas it is assumed that the in- 
sulin secretion rate depends on the intracellular Ca 2+ 
concentration. For the macroscopic description, we con- 
sider two possible insulin induced inhibitory pathways 
affecting the conductance gK(ATP) ■ In addition, electri- 
cal coupling through gap junctions between (3 cells in an 
islet is assumed. Obtained via extensive numerical simu- 
lations are a variety of oscillatory behavior of the bursting 
electrical activity, calcium concentration, and insulin se- 
cretion, which are in good agreement with those observed 
experimentally. 

There are four sections in this paper: Section II intro- 
duces the model system, described by coupled differen- 
tial equations, together with the appropriate parameters. 
The coupled equations are integrated numerically and the 
results are presented in Sec. III. Finally, a detailed dis- 
cussion, together with a brief summary, is given in Sec. 
IV. 



II. MODEL SYSTEM 

To describe the whole feedback process, we incorporate 
the bursting mechanism of f3 cells given by Sherman Q 
and the description of glucose- insulin oscillations by Maki 
and Keizer 0, • We begin with the dynamic equation 
for the time-dependent behavior of the membrane poten- 



tial Vi of the ith (3 cell 0: 

dVi 



T ~j~j~ ~ ~ ^Ca — IiK — HS ~ liK(ATP) 



Y^g c {Vi -Vj), 



(1) 



which describes current balance with the relaxation time 
t. The right-hand side of Eq. includes contributions 
from various current channels, which are all given in 
units of voltage (i.e., with the dimensionless conductance 
as below): calcium current lea, potassium current Ik, 
generic slow current Is, and background current Ik(atp) 
through the ATP-dependent potassium channel. The last 
term, where the prime in the summation stands for the 
restriction that j be a nearest neighbor of i, represents 
the electrical coupling via gap junctions of (dimension- 
less) conductance g c between nearest neighboring cells in 
an islet. 

The first three types of current depend on the mem- 
brane potential Vi and on the gating variables Nj or Si, 
and are given by 

hca = gCaMooiVi - V C a) 

I iK = gKNiiVi-Vn) 
l is = 9sSi{Vi - V K ) 

with constant (dimensionless) conductances gca, gx, and 
gs and the reversal potentials Vca and Vk for Ca 2+ and 
K + ions, respectively. The gating variables Ni and Si are 
governed by the equations 



dNi 



TS- 



dt 
dSj 
dt 



S 



(2) 
(3) 



where the activation values and Soo as well as 
in general depend on the membrane potential. Thus the 
calcium current has been assumed to respond instanta- 
neously to a change in the membrane potential via the 
voltage-dependent activation M m while the potassium 
current is governed by the gating variable Ni via Eq. J2J. 
These two currents are responsible for generating the ac- 
tion potential during the active phase of bursting. The 
generic slow current, which may corresponds to the in- 
hibitory potassium current, is gated by the slow variable 
Si via Eq. © and responsible for switching between the 
active and silent phases. We use Boltzmann-type expres- 
sions for the equilibrium values , , and of the 
voltage-dependent activation: 



Aoo — 



1 



l + exp[(Vx-Vi)/6 x ] 



where X denotes 



with appropriate constants Vx and Ox 
M, N, or S. 

On the other hand, Ik(atp), which is responsible for 
setting the plateau fraction, i.e., the ratio of the active 
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phase duration to the burst period, satisfies 



HK(ATP) 



9K(ATP)Pi{Vi - V K ), 



where the conductance Qk(atp) depends on the glucose 
concentration and pi is the opening probability of the 
K(ATP) channel. In general, an increase in the intra- 
cellular glucose concentration due to injection of glucose 
raises the ATP/ADP ratio after being metabolized; this 
in turn reduces qk(atp \ an d induces depolarization of 
the membrane potential |2Q] . Considering the time delay 
in this process, i.e., between the increase of the intracel- 
lular glucose concentration Gi n and the reduction of the 
conductance gK(ATP), we assume that gK(ATP) at time 
t is a decreasing function of Gi n at time t—t^, where the 
delay time td measures the time needed for metabolizing 
glucose transported into the (3 cell. To be specific, we 
take the Hill equation pi) : 



9k(atp) (t) 



9i - 92 



[G in (t-t d )/G K } b fD{Ji) + l 



92, 



with the appropriate Hill coefficient b [22|. Note the 
sharp contrast to the constant giaATP) m existing stud- 
ies on the bursting mechanism |l| or to the case that 
9k(atp) is taken to be a function of ATP 0. Here g\ 
and <?2 represent the maximum and minimum values of 
9k(atp), respectively, whereas Gk is the value of Gj n for 
which 9k(atp) reduces to its medium value [g\ + 52 ) / 2 . 
The function /d( J) describes possible direct effects of in- 
sulin on 9k(Atp)- Namely, while we have simply fu = 1 
in the absence of such direct effects, in their presence 
we take fo(Ji) to be a decreasing function f(Ji) of the 
inhibition variable J, : 



1 + Jo' 



1 + ( Jo/ JiY 



with the inhibition variable governed by 



tj 



dJj 



Ji. 



(4) 



The (slow) relaxation time tj and the constant Jo depend 
on the pathway, whereas is taken to be an increasing 
function of the insulin concentration H: 



Joe 



Hi 



Hj + H 

with appropriate constant Hj (see Table I). 

Although there are many sources of noise in a biological 
system, we here restrict ourselves only to the stochastic 
opening and closing of ion channels, in particular those 
of K(ATP) channels, and assume that the dynamics of 
the opening probability pi in the ith cell is described by 
the stochastic differential equation 0, |2^| 



dpi 
dt 



71 



"72 

(I- Pi) Pi +Tli{t) 



(5) 



with appropriate constants 71 and 72 and relaxation time 
t p , where i}i(t) is the Gaussian white noise with zero mean 
and variance 



7l(l - Pi) +J 2 Pi 
Tp^KiATP) 



8 lJ 5{t 



0- 



Note that as the number n K ^ TP ^ of K(ATP) channels 
per cell grows large {riK(ATP) ~ * °°), we have r\i — > 
and obtain the usual deterministic equation for pi. 

Modifying the model for a perifusion system 0, we 
now obtain equations describing glucose regulation. In 
view of the perifusion, we consider conservation of insulin 
and write the equation for the insulin concentration Hi 
around the ith cell in the form: 



dHj 
dt 



koiH - H ), 



(6) 



where R s is the rate of insulin secretion, fco is the flow 
rate for each cell, and Hq is the background concentra- 
tion. The effective volume f2 of a j3 cell is given by the 
volume of the islet divided by the number of (3 cells in 
it. Experimental observations indicate that insulin secre- 
tion from (3 cells increases with the intracellular calcium 
concentration and with the amount of insulin stored for 
rapid secretion |24) . In particular, the secretion rate is 
known to depend on the Ca 2+ concentration quartically. 
We thus assume a quartic function of the concentration 
Ci of the free intracellular Ca 2+ : 

Sr 



Rs — RoCi 

^max 

where Rq sets the scale, Ci is the rescaled (dimensionless) 
calcium concentration defined to be Ci = (C\ — Cb)/Co 
with the background value Cf, and the appropriate scale 
Co, and Sr is the storage amount of insulin prepar- 
ing rapid secretion in the readily releasable vesicle pool 
(RRVP) of the /3-cell, with the maximum value S max . 
The latter changes according to the relation 



dS H 
dt 



Rr — Rs, 



(7) 



where the refilling rate R r of insulin granules to the 
RRVP is proportional to the remaining amount: 



Sr) 



with appropriate constant a r . The Ca + concentration 
Ci in the ith cell is governed by the rate equation: 



dd 

~df = - £{aUca 



M s - M P - M N + M, 



leak ) 



(8) 



where e is the fraction of free Ca 2+ ions in the cytoso- 
lic compartments, a is the proportionality constant be- 
tween the current flow and the concentration reduction 
rate. Clearance terms, Mg, Mp, and M/v are operated 
by sarco-endoplasmic reticulum Ca 2+ -ATPase (SERCA) 
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pumps, plasma-membrane Ca 2+ ATPase, and plasma- 
membrane Na/Ca 2+ exchangers, respectively |25j, and 
given by 



M s 



M P 



]\/rmax 

M S 1 , fv.in.Y2, 



M rnax 



1 + {K s /C t f 

1 [H + ] 
1 + (Kp/d) [H+] + [KJ 



M N = fcjvCj 



with appropriate constants M™ aa; , M™ 1 , fejy, ^Cs; an d 
A'p as well as [H + ] and [K a ] adjusted to pH= 7.40 and 
pK a = 7.86 for activation of pumping by protons. The 
last term Mi ea k represents the flux out of the endoplasmic 
reticulum (ER) and is described by 

Mi ea k = £er(Cf — Cj), 

where £ er is the rate of Ca 2+ release from the ER and Cf r 
denotes the Ca 2+ concentration in the ER. The latter is 
governed by the rate equation 



dCf 
dt 



V C yt 



(M S ~ Ml ea k), 



(9) 



where e er is the fraction of free Ca 2+ ions in the ER, v cyt 
and the volumes of the cytosolic and of the ER 

compartments, respectively EH3 • 

Finally, we consider the rate equation for the glucose 
concentration Gi n in the ith cell: 



dG tr 
dt 



n + r 2 



(10) 



where r\^ 2 are the rates of glucose uptake through GLUT- 
1 and GLUT-2 transporters, respectively, and r m is the 
rate of glucose metabolism J22] . Here the difference be- 
tween the glucose concentration in blood plasma and the 
injected glucose concentration Go has been disregarded 
for convenience, with the assumption that the former sat- 
urates rather quickly to the latter. We thus take the rates 
through GLUT-1 and GLUT-2 transporters to be simple 
increasing functions of the extracellular glucose concen- 
tration Go and choose an increasing function of the in- 
tracellular glucose concentration Gi n for the metabolism 
rate: 



Cl(Go — Gin) 



Hi 



(1 + G /Kx){Kx + Gin) H r + Hi 

c 2 {Go-G m )fi(Ji) 
(l + G /K 2 )(K 2 +G m ) 

1 in 

G 

in 



K m + Gy 



The constants K%, c 2l K 2 , c m , and K m are determined 
from experimental results, whereas the constants C\ and 
H r are chosen to yield appropriate intracellular glucose 
concentrations in simulations. Note that n has been 



taken to be an increasing function of Hi, which mod- 
els the recruitment of GLUT-1 transporters by insulin 
and results in positive feedback. On the other hand, to 
take into account indirect effects of insulin on gK(ATP)i 
r 2 has been chosen to depend on the function /j of the 
inhibition variable Ji, which accounts for inhibitory ef- 
fects of insulin on its own release via the inhibition of 
GLUT-2; this assumption is made as a possible interpre- 
tation of the recent observation of the inhibiting role of 
insulin 0, |n|, [^. Accordingly, r 2 tends to decrease 
as the insulin concentration is raised. Note in partic- 
ular that unlike the form in the existing study |2j, r 2 
in the above form always remains positive, thus avoid- 
ing the unphysiological case of negative values. To com- 
plete the model, we assume the two-step process in a f3 
cell, which includes manifestation of the inhibitory sig- 
nal from insulin receptors in the presence of insulin (first 
step) and the inhibition of GLUT-2 transporters by that 
signal (second step) ; this provides the slow negative feed- 
back process, which may again be described by the de- 
creasing function /(Jj) (for the direct effects before) with 
appropriate values of the parameters. 

The whole feedback model is thus described by the 
ten coupled equations [Eqs. (JIJ - Q10[l] for the ten vari- 
ables, together with the equations giving dependence on 
those variables. It should be stressed that the bursting 
mechanism and the glucose regulation are incorporated, 
with the glucose-dependent conductance Qk(atp) and 
the calcium-dependent rate of insulin secretion as well 
as the insulin-dependent rate of glucose uptake playing 
the role of a bridge between the two. The appropriate 
values of the parameters used in the whole model are 
given in Table I. 



III. SIMULATION RESULTS 

We integrate numerically the ten equations, with the 
parameter values given by Table I, by means of the fifth- 
order Runge-Kutta method. To keep track of the state 
of every single channel is computationally expensive; it is 
sufficient to monitor only the total number of open and 
closed channels per cell. The time step At in the nu- 
merical integration should be chosen sufficiently small so 
that the probability of two channel openings (or closings) 



during one time step is negligible, say, less than 



A< < 0.1 min 



sings 
ifi: 



Jl n K(ATP) l2n K (ATP) 



In order to control efficiently the error in the integra- 
tion, we here use the adaptive step-size algorithm with 
Atmax = 1 m s, keeping the above restriction on At. This 
controls automatically the maximum of the next step size 
to achieve the predetermined accuracy. 
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TABLE I: Parameter values in the whole feedback model: The parameters without asterisks or single asterisks are fixed or 
adjustable ones, determined from experiments, respectively. Those with double asterisks are free parameters, not determined 
from experiments. 



g Ca = 3.6 
g K = 10 

*gs = 4 

a = 1.3 X 1CT 6 M/Vs 
e = 0.01 
gi = 3.0 

* M rnax = 2 46 M M/ S 

c m =4.1 x 10 _5 M/s 
ci = 5.8 x 10~ 4 M/s 



"feo = 1.8 x 10~ 2 s 
**H Q = OmM 



-l 



V Ca = 25mV 
V K = -75 mV 

*7i = 1 

*72 = 1 
e er = 0.01 

ff 2 = 0.6 

* M m aa; = 12 g ^ M / g 

*f m = 7.8 mM 
ifi = 1.4 mM 
*H r = 1.4 x 10" 2 mM 
*C = 10" 4 mM 



Vm = -20 mV 
V/v = -17mV 
*V S = -22 mV 
*r p = 0.30 s 

*Cer = 1.0 X 1Q- 4 

G K = 2.8 mM 
*k N = 84s" 1 



*7?o = 1.2 x 10 



d1/s 



C2 = 



5.3 x 10~ 4 M/s 
= 1.4 x 10" 4 mM 



*C b = 6.0 x 10 _5 mM 



6 M = 12 mV 
6 N = 5.6 mV 
*6 S = 8.0 mV 

*^if(ATP) = 10 

^ cy t = 10 fim 3 

b = 2.5 
K s = 0.27 ^M 
r = 2.0 x 10" 3 s" 
= 17 mM 
**f d = 90s 



'S, 



r = 20 ms 

A = 0.8 
*rs = 60s 
Q = 1.5 x 10" 15 m 3 

v er = 0.4 /im 3 
,nax = 3.6 x 10~ 17 mol 
K P = 0.50 /xM 
[#+] = 10" 7 ' 40 M 
[iCo] = lO" 7 ' 86 M 
**m = 4.0 




FIG. 2: Two-dimensional square lattice 



A. Indirect Pathway 

In order to consider indirect effects on gx(ATP) ( v i a 
negative feedback to GLUT2 by insulin) , we take fjj = 1 
and fi = f with Jo = 0.35, tj = 1 min, and other pa- 
rameter values in Table I. In this manner we have per- 
formed simulations of an islet, which consists of 64 cells 
arranged on a two-dimensional square lattice shown in 
Fig. |21 and obtained the membrane potential, Ca 2+ con- 
centration, and insulin secretion at the coupling conduc- 
tance g c = 0.06. Figure|21shows the obtained time depen- 
dence of the membrane potential of a single cell randomly 
chosen among 64 cells for extracellular glucose concentra- 
tion G = (a) 3, (b) 5, (c) 9, and (d) 19 mM. Observed is 
the bursting behavior of the membrane potential, forming 
periodic clusters of bursts. Also, the first burst persists 
longest, while following ones last progressively shorter, 
similar to experimental observations |24| . The total du- 
ration time of such a regularly bursting cluster grows with 
the glucose concentration, as shown in Fig. 0] Here the 
criterion for an active phase has been chosen to be the 
membrane potential higher than V — —55 mV; the dura- 
tion time T represents the total duration of such active 




FIG. 3: Bursting behavior of the membrane potential at 
several extracellular glucose concentrations: Go = (a) 3, (b) 5, 
(c)9, and (d) 19 mM. 



phases in a cluster, averaged over clusters. The intracel- 
lular glucose concentration Gi n of a single cell under the 
same conditions as in Fig. [3] is displayed in Fig. EI note 
the oscillatory behavior, with the amplitude and the pe- 
riod growing with the extracellular glucose concentration 
Go- 

Figures H3 and |H1 exhibit the corresponding be- 
haviors of the concentration of cytosolic calcium Ca 2+ , 
the insulin secretion rate, and the insulin concentration, 
respectively, under the same glucose concentrations as 
in Fig. [3] Unlike other quantities obtained for a single 
cell, the calcium concentration in Fig. [S] has been aver- 
aged over all cells. Note that oscillations of the aver- 
aged calcium concentration have periods similar to those 
of repetitive activation of the membrane potential; this 
reflects that most of the cells in an islet are well syn- 
chronized except for slight phase shift (2?|], as observed 
in simulations. It is observed that the Ca 2+ concentra- 
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t (min) t (min) 



FIG. 4: Duration time T of a bursting cluster depending on 

the extracellular glucose concentration. Also shown is the FIG. 6: Behavior of the calcium concentration at the same 
typical error bar, estimated from the standard deviation. extracellular glucose concentrations as in Fig. 3. 
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FIG. 5: Behavior of the intracellular glucose concentration at 
the same extracellular glucose concentrations as in Fig. 3. 




FIG. 7: Behavior of the insulin secretion rate at the same 
extracellular glucose concentrations as in Fig. 3. 



tion keeps increasing during the action potential firing 
and then decreases slowly after the firing stops. The in- 
sulin release occurs in the form of sharp peaks as soon as 
the Ca 2+ concentration increases in the cell (see Fig.0. 
After the rest period of about 5 min, the secretion rate 
of insulin reaches the maximum value and then dimin- 
ishes at the next burst, as the amount of insulin stored 
for rapid secretion reduces. The overall behavior of the 
insulin secretion rate agrees well with the experimental 
results Hij. 

When the flow rate is sufficiently high, on the other 
hand, the intracellular glucose concentration tends to sat- 
urate rather than to oscillate. To probe bursting pat- 
terns in this stationary state, we set the insulin flow rate 
ko = 9.2 s , maintaining previous values of other pa- 
rameters. Such a low insulin feedback condition may 
correspond to the experiment performed to obtain sta- 
tionary values of the intracellular glucose concentration 



(a) ' (b) 




t (min) t (min) 




t (min) t (min) 



FIG. 8: Behavior of the insulin concentration around a cell at 
the same extracellular glucose concentrations as in Fig. 3. 



7 




40 50 

t (min) 



40 50 

t (min) 



FIG. 9: Behavior of the membrane potential in the stationary 
state of the intracellular glucose concentration. The extracel- 
lular glucose concentrations are the same as those in Fig. 3. 
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FIG. 10: Saturated intracellular glucose concentration de- 
pending on the extracellular glucose concentration. Squares 
represent the data obtained from experiment |2Sj. whereas 
triangles are the data obtained via simulations and fitted by 
the dashed line. 



at several extracellular glucose concentrations |28] . The 
resulting behavior of the membrane potential is displayed 
in Fig. depending on the extracellular glucose concen- 
tration. The duration time of each burst is also observed 
to grow with the extracellular glucose concentration; on 
the other hand, oscillatory patterns with periods of about 
10 min do not emerge. It is thus concluded that the emer- 
gence of slow oscillations of the membrane potential de- 
pends upon the flow rate of the circulating extracellular 
fluid, i.e., blood. The lack of slow oscillations led by min- 
imal insulin feedback matches well with the experimental 
results of IRS-1 (insulin receptor substrate 1) knockout 
mice [29). 

Figure HOI shows how the intracellular glucose concen- 
tration Gin depends on the extracellular glucose concen- 
tration Go- Shown together with the experimental re- 
sults [2^ are the saturated values of Gi n obtained from 
simulations with Go varied from 1 mM to 26 mM. It is 
pleasing to see the excellent agreement between the ex- 
perimental data and the simulation data. In particular, 
the dependence of d n on Go yields reasonably good lin- 
ear fitting, with the slope 0.81. 

We also investigate the case that the concentration of 
injected glucose is varied (with k Q — 9.2 s _1 ), to probe 
time-dependent responses to sudden changes of the ex- 
tracellular glucose concentration. For simplicity, the ex- 
tracellular glucose concentration Go is given as a step 
function of time, with possible time delay disregarded. 
Figure 1111 shows the behaviors of the membrane poten- 
tial, calcium concentration, and intracellular glucose con- 
centration as glucose is injected with the concentration 
G = 5,9,12,19, and 11 mM at time t = 20,30,40,50, 
and 60 min, respectively. They are observed to display 
behaviors in close response to the change of the glucose 
concentration. 

To confirm the assertion that insulin inhibits burst- 
ing of membrane potentials, we examine the behavior of 
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FIG. 11: Behavior of (a) the membrane potential, (b) cal- 
cium concentration, and (c) intracellular glucose concentra- 
tion when glucose is injected as a step function of time: 
Go = 5, 9, 12, 19, and 11 mM at time t = 20, 30, 40, 50, and 
60 min, respectively. 
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100 nM Insulin 




FIG. 12: Behavior of the membrane potential under Go = 
10 mM, when insulin corresponding to Ho = 100 nM is in- 
jected during the time t — 10 min to 25 min. The experimen- 
tal values ci = 2.9 x 10 _4 M/s and fco = 9.2 s _1 are used. 
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the membrane potential in response to the injection of 
insulin. Assuming that the basal plasma insulin concen- 
tration changes rapidly, we represent the insulin stimulus 
by an increase in the basal insulin concentration Hq . Fig- 
ure 1 121 shows the membrane potential in response to the 
insulin injection corresponding to Hq = 100 nM during 
the time t = 10 min to 25 min under the glucose concen- 
tration Go = 10 mM. (Note that the parameter values 
for normal /3 cells are used here for comparison, differ- 
ent from those for hamster insulinoma tumor cells in Ta- 
ble I. See Ref. y and references therein.) It is observed 
that the bursting action potential disappears after about 
6 min from the beginning of the insulin injection and ap- 
pears again in about 3 min after the cease of injection. 
This is consistent with the experimental result [ill ]. 

Finally, for a more realistic model, we allow variations 
in some parameters among the cells, and examine the 
effects of such inhomogeneity. We have considered varia- 
tions in the coupling conductance g c , choosing randomly 
its value in the range 0.05 < g c < 0.07. The resulting 
behavior of the membrane potential is found to be qual- 
itatively the same as that for constant g c . 



FIG. 13: Bursting behavior of the membrane potential via 
the direct pathway mechanism with tj = 3 min, at the same 
extracellular glucose concentrations as in Fig. 3. 
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FIG. 14: Behavior of the intracellular glucose concentration 
at Go = 9mM and tj = 3 min in the presence of the direct 
pathway. The inset shows an enlarged view. 



B. Direct Pathway 

There are reports that insulin itself inhibits insulin se- 
cretion in /3-cells although the exact mechanism is not 
clear [jjj OH. Specifically, experimental data seem to 
suggest that the negative feedback of insulin acts directly 
on ATP-dependent K + channels, suppressing action po- 
tential firing [llj . To examine this direct pathway, we set 
fj = 1, f D = f with J = 0.15, and Hj = 1.0 x 10^ 3 mM, 
varying the relaxation time tj. Figure lT3l exhibits the re- 
sulting oscillatory behaviors of the membrane potential 
for tj = 3 min, depending on the extracellular glucose 
concentration Go- The average duration of bursting is 
rather insensitive to Go, approximately 16 seconds for a 
wide range of Go- Compared with the results of the in- 
direct pathway, the duration time is short and the dose- 
response to glucose is also weak in the direct pathway. 
On the other hand, the average period between bursts 
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FIG. 15: Bursting period t p depending on tj at Go = 10 mM. 
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FIG. 16: Duration time T of a bursting cluster versus the time 
delay Td, with typical error bars estimated from the standard 
deviation. The insets exhibit the bursting patterns at Td = 
1.2 min (upper left) and 9.0 min (lower right). 



decreases with Go, which results in more insulin secre- 
tion at higher glucose levels. Unlike the indirect path- 
way, however, the clustering of bursts does not emerge in 
the direct pathway. Since the secreted insulin affects the 
ATP-dependent K + channels directly and rapidly, there 
is not enough time to activate the mechanism for burst- 
ing clusters. The behavior of the intracellular glucose 
concentration Gt n is shown in Fig.^J for the extracellu- 
lar glucose concentration Go = 9 mM and the relaxation 
time rj = 3 min. No prominent slow oscillations are 
observed (compare with Fig. |5j) . Figure 1151 shows how 
the average bursting period t p depends on the relaxation 
time rj: Observed is a rapid increase of t p , followed by 
saturation to 2 min, as tj is raised. 



To obtain bursting clusters via the direct pathway, one 
needs another formal time delay and accordingly, an ad- 
ditional process in insulin inhibition. As an attempt, 
we consider delay Td in the direct activation factor, i.e., 
foit — Td) on the conductance of the ATP-dependent K + 
channel, and show the results in Fig. It is observed 
that the total duration time T of a bursting cluster grows 
almost linearly with the delay Td- When Td is small, we 
have T ~ and there is no regular bursting patterns. 
When Td is raised to 1 .2 min, clustering patterns of bursts 
begin to appear; the corresponding behavior of the mem- 
brane potential is displayed in the upper inset panel. As 
Td is increased further, the bursting duration time keeps 
growing and the bursting pattern at Td = 9.0 min be- 
comes similar to that via the indirect pathway shown 
in Fig. OJc). Nevertheless it is rather difficult to recog- 
nize this direct pathway as an appropriate mechanism for 
bursting clusters since such large delayis apparently not 
consistent with experimental results 



IV. DISCUSSION 

Our model has employed the glucose-insulin feedback 
for slow dynamics, responsible for slow oscillations in the 
system. Unlike the fast oscillations of bursting action 
potentials, however, the origin of the slow oscillations 
with periods about 5 to 10 minutes is not understood 
well and still in controversy. One of the main difficulties 
lies in identifying the primary oscillation among many 
secondary ones. Indeed in a pancreas, many chemical 
ingredients, e.g., Ca 2+ , O2, ATP, glucose, and insulin, 
show oscillator y b ehaviors with periods of several min- 
utes 0, E3) UH> HOI • Recently, it has been reported that 
glycolysis produces spontaneous oscillations in /3-cells: 
Here Ca 2+ acts as a mediator, transducing oscillatory 
metabolism into oscillatory secretion [3(j ■ After comple- 
tion of the first version of the manuscript, we became 
aware of the recent work, which addresses the multi- 
ple bursting mode in pancreatic islets, with glycolysis 
adopted as the slow dynamics |33| • In the present study, 
on the other hand, we have proposed another possible 
mechanism for slow oscillations, namely, the negative 
feedback through the product (i.e., insulin) in response 
to the stimulus (glucose). 

In combining the slow dynamics (glucose-insulin feed- 
back) with the fast dynamics (bursting action poten- 
tials), we have allowed both dynamics to interact bidi- 
rectionally. As a result, there emerge mixed oscillations 
of both time scales, clusters of bursts, when the two dy- 
namics complement each other via the indirect pathway. 
This also keeps parallel with the glycolytic model |31| . 
and contrasts with the unidirectionally interacting model 
where slow dynamics of glycolysis influences fast dynam- 
ics one-sidedly [32|. 

We now discuss the plausibility of our model, in com- 
parison with others. Insulin secretion is synchronized 
among the (3 cells in an islet, which has been observed 
in the calcium imaging experiment p7^ . In addition, 
the pulsatile insulin secretion in a group of islets indi- 
cates additional synchronization among islets in a pan- 
creas; otherwise the secreted insulin from many islets, 
having different phases from each other, would produce 
essentially a constant level of insulin |3j| , although there 
may exist a source of some controversy [34]]. Such syn- 
chronization among islets may not be explained by the 
model adopting glycolysis as the slow dynamics, where no 
communication is available between islets. An expanded 
model was also proposed, incorporating glucose release 
in the liver: Secreted insulin causes oscillations of the 
plasma glucose through the liver, which in turn entrain 
all islets with its characteristic frequency |33j. This is, 
however, in disagreement with pulsatile insulin secretion 
at constant gluc ose concentration observed in vitro |35| 
and in vivo [36j. On the other hand, in our feedback 
model, secreted insulin from an islet affects other islets 
(via the paracrine interaction) as well as the islet itself 
(via the autocrine interaction), thus inherently allowing 
synchronization among islets even at a constant plasma 
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glucose level. Accordingly, our model can explain the 
sustained insulin oscillations at consta nt g lucose concen- 
tration observed in vitro and in vivo |35i l36j . Further, 
the recent observation that defects in IRS-1, associated 
with the reduced expression of the SERCA protein, re- 
sult in lack of slow calcium oscillations and reduction of 
insulin secretion |29j is also supportive of the insulin feed- 
back mechanism. Note, however, that the molecules for 
negative feedback are not necessarily limited to insulin 
and other products, e.g., GAB A, can also provide the 
feedback pathway |33 . 

In constructing a model for such complex biological 
systems, one may not expect perfection and impeccable- 
ness, and there are certainly some drawbacks in our feed- 
back model as well. First, our model predicts that both 
the amplitude and the period of the insulin oscillation 
increase with the plasma glucose concentration, which 
apparently does not agree with the rep ort that mainly 
the amplitude tends to change |18|, |36( . Secondly, slow 
oscillations are still observed in experiment on single j3 
cells, even when the flow rate of the local perfusion is too 
fast for secreted insulin to reach the insulin receptor and 
to operate the feedback mechanism j^l ■ In this case the 
other model adopting glycolysis as slow dynamics is still 
applicable since the glycolytic oscillation may still occur 
without the feedback through a product like insulin |3l| . 

In summary, we have incorporated the macroscopic de- 
scription of glucose regulation and the microscopic mech- 
anism of bursting behavior of (3 cells, to establish a 
model that displays inherently the observed oscillations 
of the membrane potential, cytosolic calcium concentra- 
tion, and insulin secretion in pancreatic islets. In view of 
the experimental observations, the ATP-dependent con- 
ductance gK(ATP) °f the K + channel has been taken as a 
decreasing function of the glucose concentration whereas 
the insulin secretion rate has been given by an increasing 
function of the intracellular Ca 2+ concentration. Con- 
sidered here are two possible insulin induced inhibitory 



pathways, affecting the conductance gK(ATP) ■ Further, 
f3 cells in an islet have been considered to be coupled 
electrically with their nearest neighbors with appropri- 
ate coupling conductance. By means of extensive nu- 
merical simulations, we have obtained bursting electrical 
activities, calcium concentration, and insulin secretion, 
which are consistent with those observed experimentally. 
In addition to the stationary-state behaviors, we have 
also examined how the intracellular glucose concentra- 
tion changes with the extracellular glucose concentration, 
and explored the behaviors of the system as the injected 
glucose concentration is varied successively. Then the 
inhibition effects of insulin on the bursting action po- 
tential have been probed. Finally, we have compared 
two possible pathways to the feedback modulations of 
9k(Atp) by secreted insulin. With appropriate parame- 
ters, the indirect feedback mechanism has been found to 
generate bursting clusters of action potentials. On the 
other hand, the direct mechanism in general generates 
only single bursts; it requires, e.g., rather unrealistically 
long delay time to produce bursting clusters. 

Our model may thus serve as a useful framework to pre- 
dict insulin secretion modulated by the molecular mecha- 
nism in p cells and the corresponding glucose regulation. 
In particular, our model, based on the glucose-insulin 
feedback, can be a good complement to the glycolytic 
model. Extensions toward a more realistic description of 
the biological situation are left for further study. 
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